#!/usr/bin/env python
"""

    mapping_quality_read_depth_of_uncallable.py
    [--log_file PATH]
    [--verbose]

"""

################################################################################
#
#   mapping_quality_read_depth_of_uncallable.py
#
#
#   Copyright (c) 11/3/2010 Leo Goodstadt
#
#   Permission is hereby granted, free of charge, to any person obtaining a copy
#   of this software and associated documentation files (the "Software"), to deal
#   in the Software without restriction, including without limitation the rights
#   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#   copies of the Software, and to permit persons to whom the Software is
#   furnished to do so, subject to the following conditions:
#
#   The above copyright notice and this permission notice shall be included in
#   all copies or substantial portions of the Software.
#
#   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
#   THE SOFTWARE.
#################################################################################

import sys, os

# add self to search path for testing
if __name__ == '__main__':
    exe_path = os.path.split(os.path.abspath(sys.argv[0]))[0]
    module_name = os.path.split(sys.argv[0])[1]
    module_name = os.path.splitext(module_name)[0];
else:
    module_name = __name__

# Use import path from <<../python_modules>>
if __name__ == '__main__':
    sys.path.append(os.path.abspath(os.path.join(exe_path,"/home/lg/python_modules")))



#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   options


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


if __name__ == '__main__':
    from optparse import OptionParser
    import StringIO

    parser = OptionParser(version="%prog 1.0", usage = "\n\n    %prog [options]")
    parser.add_option("-i", "--input_path", dest="input_path",
                      metavar="PATH",
                      type="string",
                      help="Path of input files. "
                          "Defaults to reading from STDIN.")

    #
    #   general options: verbosity / logging
    #
    parser.add_option("-v", "--verbose", dest = "verbose",
                      action="count", default=0,
                      help="Print more verbose messages for each additional verbose level.")
    parser.add_option("-L", "--log_file", dest="log_file",
                      metavar="FILE",
                      type="string",
                      help="Name and path of log file")
    parser.add_option("-R", "--R_plot_cmds_file", dest="R_plot_cmds_file",
                      metavar="FILE",
                      type="string",
                      help="Name and path of file to hold commands for cummulative plot by mapping quality in R")
    parser.add_option("--skip_parameter_logging", dest="skip_parameter_logging",
                        action="store_true", default=False,
                        help="Do not print program parameters to log.")
    parser.add_option("--debug", dest="debug",
                        action="count", default=0,
                        help="Set default program parameters in debugging mode.")
    parser.add_option("--mouse_euchromatin_size", dest="mouse_euchromatin_size",
                      metavar="INTEGER",
                      type="int",
                      default = 2558525933,
                      help="Size of mouse euchromatic genome 2,558,525,933bp in mm9")




    # get help string
    f =StringIO.StringIO()
    parser.print_help(f)
    helpstr = f.getvalue()
    original_args = " ".join(sys.argv)
    (options, remaining_args) = parser.parse_args()

    if not options.log_file:
        options.log_file                = module_name + ".log"

    #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
    #                                             #
    #   Debug: Change these                       #
    #                                             #
    #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    if options.debug:
        options.verbose                 = 5
        if not options.input_path:
            options.input_path              = "/data/mus/visitors/lg/uncallable/*.nc"
            original_args += " --input_path %s" % options.input_path
        if not options.R_plot_cmds_file:
            options.R_plot_cmds_file            = module_name + ".R"
            original_args += " --R_plot_cmds_file %s" % options.R_plot_cmds_file
    #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
    #                                             #
    #   Debug: Change these                       #
    #                                             #
    #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

    #
    #   mandatory options
    #
    mandatory_options = ["input_path"]
    def check_mandatory_options (options, mandatory_options, helpstr):
        """
        Check if specified mandatory options have been defined
        """
        missing_options = []
        for o in mandatory_options:
            if not getattr(options, o):
                missing_options.append("--" + o)

        if not len(missing_options):
            return

        raise Exception("Missing mandatory parameter%s: %s.\n\n%s\n\n" %
                        ("s" if len(missing_options) > 1 else "",
                         ", ".join(missing_options),
                         helpstr))
    check_mandatory_options (options, mandatory_options, helpstr)


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   imports


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#from json import dumps
from collections import defaultdict
from glob import glob
from commify import commify



#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Functions


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Logger


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

if __name__ == '__main__':
    import logging
    import logging.handlers

    MESSAGE = 15
    logging.addLevelName(MESSAGE, "MESSAGE")

    def setup_std_logging (logger, log_file, verbose):
        """
        set up logging using programme options
        """
        class debug_filter(logging.Filter):
            """
            Ignore INFO messages
            """
            def filter(self, record):
                return logging.INFO != record.levelno

        class NullHandler(logging.Handler):
            """
            for when there is no logging
            """
            def emit(self, record):
                pass

        # We are interesting in all messages
        logger.setLevel(logging.DEBUG)
        has_handler = False

        # log to file if that is specified
        if log_file:
            handler = logging.FileHandler(log_file)
            handler.setFormatter(logging.Formatter("%(asctime)s - %(name)s - %(levelname)6s - %(message)s"))
            handler.setLevel(MESSAGE)
            logger.addHandler(handler)
            has_handler = True

        # log to stderr if verbose
        if verbose:
            stderrhandler = logging.StreamHandler(sys.stderr)
            stderrhandler.setFormatter(logging.Formatter("    %(message)s"))
            stderrhandler.setLevel(logging.DEBUG)
            if log_file:
                stderrhandler.addFilter(debug_filter())
            logger.addHandler(stderrhandler)
            has_handler = True

        # no logging
        if not has_handler:
            logger.addHandler(NullHandler())


    #
    #   set up log
    #
    logger = logging.getLogger(module_name)
    setup_std_logging(logger, options.log_file, options.verbose)


    #
    #   log programme parameters
    #
    if not options.skip_parameter_logging:
        programme_name = os.path.split(sys.argv[0])[1]
        logger.info("%s %s" % (programme_name, original_args))

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Main logic


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
if __name__ == '__main__':

    total_bp_by_mapping_threshold = defaultdict(int)
    total_bp_beyond_read_depth    = 0

    bp_beyond_mapping_threshold = defaultdict(int)
    bp_beyond_read_depth      = defaultdict(int)
    bp_in_borderline_regions = defaultdict(int)

    cnt_strains = 0
    sys.stderr.write("\n\n")
    logger.log(MESSAGE, "Files from %s", options.input_path)
    for input_file_path in glob(options.input_path):

        strain_name = os.path.splitext(os.path.basename(input_file_path))[0]
        cnt_strains += 1

        cnt_lines = 0

        # progress
        sys.stderr.write("\t[")

        input_file = open(input_file_path)

        bp_per_mapping_quality = defaultdict(int)
        bp_per_read_depth      = defaultdict(int)

        for line in input_file:
            cnt_lines += 1
            if cnt_lines % 100000 == 0:
                sys.stderr.write(".")
            if not len(line):
                continue
            if line[0] == '#':
                continue
            data  = line[:-1].split("\t")
            if len(data) != 5:
                print >>sys.stderr, "WARNING: %d fields on line#%d" % (len(data), cnt_lines)
                continue
            contig = data[0]
            beg, end, mapping_quality, read_depth = map(int, data[1:])
            #read_depth > 200 mq < 40:
            if read_depth > 150:
                bp_per_read_depth[read_depth] += end - beg + 1
            elif mapping_quality < 40:
                bp_per_mapping_quality[mapping_quality] += end - beg + 1
            else:
                bp_in_borderline_regions[strain_name] += end - beg + 1
                #print >>sys.stderr, "Why are we here :%s" % line[:-1]

        bp_beyond_mapping_threshold[strain_name] = sum(bp_per_mapping_quality.values())
        bp_beyond_read_depth       [strain_name] = sum(bp_per_read_depth.values())

        for k,v in bp_per_mapping_quality.iteritems():
            total_bp_by_mapping_threshold[k] += v
        total_bp_beyond_read_depth    += bp_beyond_read_depth[strain_name]

        # progress
        sys.stderr.write("]\n")

        #
        #   Histogram per input files (e.g. per strain)
        #
        #for k in sorted(bp_per_mapping_quality.keys()):
        #    print "mq=%6d, %8d" % (k, bp_per_mapping_quality[k])
        #
        #for k in sorted(bp_per_read_depth.keys()):
        #    print "rd=%6d, %8d" % (k, bp_per_read_depth[k])
        #


    #
    #   Results summary per input files (e.g. per strain)
    #
    logger.log(MESSAGE, "%d strains" % cnt_strains)
    logger.log(MESSAGE, "{0:>25s} {1:^25s} {2:^25s} {3:^25s}".format("Strain name", "Total", "Mapping threshold", "Read depth threshold", "Borderline"))
    logger.log(MESSAGE, "{0:>25s} {1:^25s} {2:^25s} {3:^25s}".format("=" * 20, "=" * 25, "=" * 25,"=" * 25,"=" * 25))
    inv_mouse_euchromatin_size = 100.0 / options.mouse_euchromatin_size
    for strain_name in sorted(bp_beyond_mapping_threshold.keys()):
        total_bp = sum((bp_beyond_mapping_threshold[strain_name], bp_in_borderline_regions[strain_name], bp_beyond_read_depth[strain_name]))
        logger.log(MESSAGE, "%25s %13s %8.1f%% %13s %8.1f%% %13s %8.1f%% %13s %8.1f%%" %
                            (   strain_name                             ,
                                commify(total_bp), total_bp* inv_mouse_euchromatin_size,
                                commify(bp_beyond_mapping_threshold[strain_name]), bp_beyond_mapping_threshold[strain_name]* inv_mouse_euchromatin_size,
                                commify(bp_beyond_read_depth[strain_name]       ), bp_beyond_read_depth[strain_name]       * inv_mouse_euchromatin_size,
                                commify(bp_in_borderline_regions[strain_name]   ), bp_in_borderline_regions[strain_name]   * inv_mouse_euchromatin_size,
                            ))


    #
    #   Cummulative proportion by mapping quality averaged across all input files (e.g. all strains)
    #
    inv_mouse_euchromatin_size = 100.0 / options.mouse_euchromatin_size / cnt_strains
    logger.log(MESSAGE, "")
    logger.log(MESSAGE, "Average %.2f%% above read depth threshold" % (total_bp_beyond_read_depth * inv_mouse_euchromatin_size))
    running_total = 0
    logger.log(MESSAGE, "")
    logger.log(MESSAGE, "Proportion below mapping threshold on average")
    logger.log(MESSAGE, "%20s %10s" % ("Mapping Quality", "Cummulative Proportion"))
    logger.log(MESSAGE, "%20s %10s" % ("=" * 18, "=" * 22))
    proportions = []
    for k,v in sorted(total_bp_by_mapping_threshold.iteritems()):
        running_total += v
        proportion = running_total * inv_mouse_euchromatin_size
        proportions.append(proportion)
        logger.log(MESSAGE, "%20d %10.2f%%" % (k,proportion))


    r_format_proportions = "c(%s)" % (",".join(map(str, proportions)))

    R_file = open(options.R_plot_cmds_file, "w")
    R_file.write("""
    # create jpeg file
    jpeg(file="uncallable_by_mapping_quality.jpg", width = 1200, height = 1200, units = "px",)
    uncallable = 0.01 * %s
    mq = 0:39
    plot(mq,uncallable, type='b', xlab = 'Mapping quality', ylab = "Cummulative uncallable proportion", main = "Cummulative proportion of uncallable genome by mapping quality" )
    # plot onto axes
    points(c(-50, mq[11],  mq[11]), c(uncallable[11], uncallable[11], -1), type = 'l', col = "blue")
    points(c(-50, mq[21],  mq[21]), c(uncallable[21], uncallable[21], -1), type = 'l', col = "blue")
    points(c(-50, mq[31],  mq[31]), c(uncallable[31], uncallable[31], -1), type = 'l', col = "blue")
    dev.off()

    # create svg file
    svg(file="uncallable_by_mapping_quality.svg")
    uncallable = 0.01 * %s
    mq = 0:39
    plot(mq,uncallable, type='b', xlab = 'Mapping quality', ylab = "Cummulative uncallable proportion", main = "Cummulative proportion of uncallable genome by mapping quality" )
    # plot onto axes
    points(c(-50, mq[11],  mq[11]), c(uncallable[11], uncallable[11], -1), type = 'l', col = "blue")
    points(c(-50, mq[21],  mq[21]), c(uncallable[21], uncallable[21], -1), type = 'l', col = "blue")
    points(c(-50, mq[31],  mq[31]), c(uncallable[31], uncallable[31], -1), type = 'l', col = "blue")
    dev.off()
                 """ % (r_format_proportions,r_format_proportions))



